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PREJTJD - A COMPUTER CODE FOR THE PRELIMINARY ANALYSIS OF 


TWO-DIMENSIONAL PULSE HEIGHT ANALYZER DATA 
by Thor T. Semler 
Lewis Research Center 

SUMMARY 

A computer program, PREJUD, has been developed for the processing of two- 
dimensional analyzer, pulse height data. The program provides an automated separation 
of gamma counts from neutron counts and produces data in a format suitable for input into 
spectrum unfolding codes, for example, FERDOR. 

A strategy for the separation of gamma and neutron counts is detailed. Data smooth- 
ing techniques and a spline interpolation technique are described in some detail. 

An example of the unfolding of data from a plutonium -beryllium source is shown to 
illustrate the use of the code. 


INTRODUCTION 

It has been found experimentally that the overall shapes of the scintillation pulse in- 
duced in trans-stilbene and other organic scintillators by gamma-rays and fast neutrons 
differ (ref. 1). This effect has been used as the basis of a technique (pulse shape dis- 
crimination) demonstrated by Brooks (ref. 2) to use organic scintillators to discriminate 
between gamma-ray and neutron induced scintillations. Details of the circuitry developed 
for pulse shape discrimination can be found in reference 1. 

At the Lewis Research Center a modified Owen pulse shape discrimination circuit 
has been used (ref. 3). The output of the pulse shape circuit and a linear output from the 
photomultiplier tube are amplified and recorded in a 4096 channel two-dimensional ana- 
lyzer. The proton recoil and Compton recoil spectra are measured and the FERDOR 
code (ref. 4) is used to unfold the proton recoil spectra. 

Generally, the preparation of data for the FERDOR code has been done manually 
(ref. 5). The PREJUD ( pre liminary adjustment of data) code described herein allows 
one to prepare suitable input to FERDOR, or other proton recoil spectra unfolding codes, 



from the multichannel analyzer data with up to four different gain settings and a minimum 
of inspection. 


PROBLEM ANALYSIS 

The data from a single experiment, one gain setting and one counting time, are pro- 
vided as a two-dimensional array, 64 by 64, of count rate as a function of energy. The 
results of a well resolved experiment for one gain setting are shown in figure 1. The 
first phase of the data analysis is the separation of the gamma induced counts from the 
neutron induced counts to provide a set of neutron counts as a function of energy. The 
technique used is described in the next section, MINIMUM FINDING STRATEGY. Next 
the process is repeated for the background, if any, which is then subtracted from the 
data. If more than one gain setting has been used, the previous steps are repeated to 
provide several sets of neutron induced counts as a function of energy. Next these sets 
of data are smoothed to remove some of the experimental scatter. The smoothing tech- 
niques are derived in the section SMOOTHING TECHNIQUES. The sets of smoothed data 
are then normalized to the same gain setting and merged. The merged data are smooth- 
ed, and then these results are interpolated and extrapolated at specified energy points so 
as to provide input to the code FERDOR (ref. 4). The interpolation is described in the 
section SPLINE INTERPOLATION. The extrapolation techniques (low energy and high 
energy) are given in the section EXTRAPOLATION OF DATA. The error estimate given 
these interpolated and extrapolated values is detailed in the ERROR ANALYSIS section. 
An example of the use of the codes PREJUD and FREDOR is shown in the RESULTS sec- 
tion. The data input for the program is given in appendix A and the program listing is 
given in appendix B. 


MINIMUM FINDING STRATEGY 

As indicated in the section on PROBLEM ANALYSIS, the first task of the code is the 
separation of counts into electron recoil counts caused by gamma rays and proton recoil 
counts caused by neutrons. This is accomplished by searching along rows of equal energy 
for the first minimum and summing to the left for electron recoil counts and to the right 
for proton recoils. Half of the minimum is added to each sum. An idealized situation 
is shown in figure 2, where the minimum between the Compton recoil spectra and the 
proton recoil spectra is easily discerned. The technique used to separate the two sets of 
counts is described in this section and shown in flow chart form in figure 3. 

(1) A row of counts A. (i = I , 64) corresponding to the same energy is read into the 
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minimum finding subroutine, where I is the number of columns ignored, if any. 

(2) If a minimum i m - n has been found in an earlier row of this data set, start the 
minimum search five positions to the left of the earlier minimum and then go to step (12). 

(3) If the magnitude of any of the first four data cells is greater than 100, go to 
step (11). 

(4) Check for zeros in the first three data cells. If all are zero, go to step (9). 

(5) Set i to two larger than the first data cell in this row. Check for two sequentially 
zero data cells. If they are found, go to step (7). 

(6) Increment i counter by one and return to step (5). 

(7) Subtract 1 from i and sum from left to i for electron recoil counts and from i 
to right for proton recoil counts. 

(8) Set i m ^ n equal to i and return to calling routine. 

(9) Sum entire row since there are no electron recoil counts and set equal to proton 
recoil counts. 

(10) Set i fflin equal to first data cell i and return to calling routine. 

(11) Set i equal to one larger than first data cell i. 

(12) Check for condition A i-1 > A i < A i+1 : if this condition is satisfied, go to 
step (17). 

(13) Check for condition A i- 2 > A i-i = A i< A i + r if this condition is satisfied, go to 
step (17). 

(14) Check for condition A i -3 ^ A i-2 A i-1 ~ A i < A i+1 ; if this condition is satisfied, 
go to step (17). 

(15) Check for A i = 0; if it is true, set i min = i and go to step (19). 

(16) Increment i by 1 and go to step (12). 

(17) If A i > (A. +1 - 3 y/A i+1 ), go to step (22). 

(18) If A i > A i+2 , go to step (16). 

(19) Set i - n = i and perform double summing. 

(20) If the sum of the proton recoil counts is less than seven times the counts in 

A i min or ^ the ra ^° °f electron recoil counts over proton recoil counts is greater than 

10^, the sums and i . are set equal to zero. 

min n 

(21) Return to the calling routine. 

(22) If any of the differences (A i+4 - A^), (A i+3 - A^, or (A i+2 - A^ is negative, a 
local minimum has not been found and the code returns to step (16). 


SMOOTHING TECHNIQUES 

Once the data have been resolved and stored in bins corresponding to the channel 

energy, they may be smoothed prior to normalization and spline fitting. Comment cards 
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are inserted in the code listing (appendix B) indicating how one may either smooth the 
data more or may eliminate most of the smoothing. A five point least squares fit of the 
log counts has been chosen to smooth these data. The log counts have been chosen so as 
not to accentuate large count rates in the smoothing. Since the curvature may well be 
changing considerably between data points, a five point least squares parabola has been 
constructed to smooth the equal interval data. 

Begin by taking five points centered on (x Q , y Q ), that is, (x_ 2 , y_ 2 ), (x_j, y_j), (x Q , 
Yq) , (xj, yj), and (x 2 , y 2 ), and transform the abscissa values so that the transformed 
values are (-2, y_ 2 ), (-1, y_ 1 )(0, y Q ), (1, y x ) , and (2, y 2 ). What then is y(x Q ), the 
smoothed least squares value of yg? Let the approximating parabola be given by equa- 
tion (1): 

y = 3-g + a^Xrp + a 2 X T (1) 

where X^, is the transformed X value. 

Thus, the problem has been reduced to minimizing S, with respect to the coefficients 
ag, ap and a 2> where S is given by equation (2) 

S = Y ( y i " a 0 ■ a l X Ti " a 2 X Ti) 
i=2 ' ' 

Hence, the normal equations are 
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Since the following relations hold, 
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i=~2 



i=-2 





i=-2 


34 


we may simplify the normal equation in the following manner: 


|r- = £ ^ " 5a o " 10a 2 = 0 

5a 0 i=-2 


as 

9a, j 


2 

' £ - IOa i 


i=-2 


= 0 


_as 

3a 


= £ (4in) - 10a o - 34a 2 = 0 


2 i=-2 


The value of the parabolic equation (1) obtained at X T = 0 is ag; therefore, one 
may solve the system of preceding equations for ag: 

-3.4[(y_ 2 + y_ 1 + y Q + y x + y 2 ) - 5a Q - 10aJ = o 
(4y„ 2 + y_i + Yi + 4y 2 ) - 10a Q - 34a 2 = 0 


Adding the two previous equations and solving for ag gives the following expression: 
y(at x Q ) = a 0 = A- ( -3y _ 2 + 12y_ 1 + 17y Q + 12yj - 3y 2 ) 


A smoothing of this variety with endpoint corrections (ref. 6) is performed on the 
log counts prior to the normalization of the data and the spline fitting of the results. 

Once the results of measurements for several gain settings have been normalized to 
the same gain setting and merged into the same data set, a least squares weighted ex- 
ponential three point smoothing is performed on the merged data set (ref. 7). The mid- 
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point is weighted by unity and the two surrounding points are weighted by the ratio of 
their errors to the central point error. This final smoothing produces a smooth set of 
data suitable for spline fitting and interpolation. 


SPLINE INTERPOLATION 


Once the data points have been obtained in a suitable form, normalized to same gain 
setting, and merged, the data must be interpolated at arbitrary energy points to provide 
input data to a spectrum unfolding code, such as FERDOR. A spline interpolation has 
been chosen for this code, so as not to introduce extraneous undulations into the proton 
recoil spectrum, as a least squares polynomial fitting of the data might. 

A spline is a device used in drafting for drawing smooth curves; it consists of a 
flexible strip of wood or metal which can be bent to pass through any points. If the spline 
is uniformly flexible and if it lies along the plane curve y = F(x), its potential energy is 


f [f"(x)] 2 [l + F'(x) 2 ] 


-3/2 


dx 


( 3 ) 


where A is a proportionality constant. The form of the spline is such as to minimize 

equation (3) consistent with the constraints imposed by the points it passes through. A 

2 

mathematical spline is defined in the same manner except that the F’(x) term is re- 
moved. 

A cubic spline has the following three features: 

(1) If there exists a set of numbers a = Xq < . . . . < = b and a set of cor- 
responding y values y 0 , . . . , y N , the cubic spline F(x) satisfies F(X K ) = y R 

(K= 0, 1 N). 

(2) F(X) and F'(X) are continuous (Xq, X n ). 

f X N r 12 

(3) J F"(X) dx exists in the interval (Xq, Xj^) and is minimized subject to 

Xq 

the two features previously given. The cubic spline is constructed in the following manner 

(refs. 8 and 9): Let S. = X., -X. and t. = F"(X.) in the interval (X., X. ,). F(X) is 

3 J +i J 3 3 3 3 +i 

then 
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F(X 


[(X - x j)y itl + (X H 1 - X) yj 


s i 


(S- + x. , -x)t. + (S. + X - X.)t iH 
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] 3+1 6S. 

J 


The condition that F'(X) is continuous at Xj is expressed as 


F’(X.) 


(y i+i - y i> . (2t i * Vi )s j . < y j - y j-i> ^ < 2t j - 


s, 


s, 


i-i 


0 = 2, 3 N - 1) t N = tj =0.5 


tion 


The values of tg, . . . , t n _j are determined from the solution of the matrix equa- 


where 


At = u 


'2(Sj + S 2 ) 
S 


o 


o 


2 2(S 2 + S3) S 3 


X 
X X 


X 


X 


\ \ x 

^ 2 -^N-2 4 S N-1* 


th 

and the j n element of u is given by 


u i = 6(y i+ i 


y i )/s r 6 (y j- y i-i )/s < 


j-1 


The coding details are given in reference 9. 

The spline interpolation thus provides a smooth interpolating function, without the 

difficulties introduced by polynomial curve fitting; for example, polynomial curve fitting 
often introduces extraneous detail into the data, undulations not found in the data (ref. 10). 
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EXTRAPOLATION OF DATA 


The values at which interpolated data are required are fixed in the FERDOR input. 

At times the experimental values do not span the energy region required by FERDOR. 
Thus, the code PREJUD has provision for both low-energy and high-energy extrapolation. 

In both regions (low and high energy) weighted least squares exponential curve fits have 
been used (ref. 7). 

For the low-energy extrapolation the first three low-energy proton recoil counts are 
weighted 10. 0, 1.0, and 0. 1 and an exponential is fitted in the weighted least squares 
manner and the FERDOR low-energy points are thus extrapolated. 

The high-energy extrapolation uses the last five points with more than 25 counts per 
point and weights them 1, 2, 4, 8, and 16. One may easily change the value to a still 
smaller number of counts or larger as noted in the comment cards. Again a weighted 
least squares exponential fit is performed and used to extrapolate the data. 

The extrapolations join the spline approximation in a smooth manner because of the 
rather large weights given end points. 


ERROR ANALYSIS 

As the FERDOR code requires not only counts, at previously decided upon energies, 
but count errors as well, the code PREJUD computes the error associated with each 
energy point prior to normalization. It computes the square root of the number of counts 
as the standard deviation of each point. If the data are normalized, the errors are nor- 
malized as well. 

In order to ascertain the error of one of the interpolated points, the code determines 
the errors of the two surrounding points and uses a linear interpolation of the errors. 

The errors associated with the low-energy extrapolation are in the same ratio as the 
lowest energy data point. That is, if the lowest energy data point is 10 000 counts with 
an error of 100 counts, an extrapolated point of 50 000 counts would have an error of 
500 counts. 

The errors associated with the high-energy extrapolation are the same as the error 
of the last point used in the extrapolation. That is, if the last point used is 36 counts with 
an error of 6 counts, an extrapolated point of 3 counts will be assigned an error of 6. 


RESULTS 

Example of Use of PREJUD 

A Pu-Be source has been measured with a 1- by 1-inch (2. 54- by 2. 54-cm) liquid 
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organic scintillator operated at three gain settings, with correction for background. A 
separation of the proton recoil spectrum has been performed manually. Furthermore the 
input for FERDGR has been interpolated by hand. This proton-recoil spectrum has been 
unfolded by FERDOR to obtain a neutron spectrum. The same data have been treated by 
the code PREJUD and unfolded by FERDOR. The results of the manual -FERDOR and 
PREJUD-FERDOR are compared in figure 3. The unfolded neutron spectra from these 
two analyses are shown as two cross-hatched regions one standard deviation about their 
respective mean spectra. 


Discussion of Example Problem 

It is found that the two spectra overlap (are statistically equivalent) over the range 
of neutron energies from about 1 to 10 MeV. Given more analyzer channels in the low- 
energy region and adequate counting statistics in the high-energy region, the results of 
hand analyses and PREJUD-FERDOR unfoldings should match over the entire energy 
range. 

The reason for the rather spread of the PREJUD-FERDOR spectrum is the rather 
conservative estimate of error; that is, PREJUD uses a linear interpolation of error 
values. The error values estimated in hand analyses are apparently associated with the 
smaller error of the two surrounding measured points. 

The hand calculation represents on the order of at least 1 man day of work; whereas 
the execution time of PREJUD for this case was 0. 54 minute on the Lewis IBM 7094-11. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, October 21, 1970, 

129-02. 
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APPENDIX A 


INPUT FORMAT FOR PREJUD 

Card 1 Card columns 2 to 80 - Identification any valid Hollerith symbols 

Card 2 Card columns 1 to 5 - format 15, number of different gain settings fixed point 

number from 1 to as many as 4 

Card 3 Card columns 1 to 30 - format (2F10. 5, 15, 15), gain setting information per- 
taining to following set of data: 

Columns 1 to 10 (F10. 5) Columns 11 to 20 (F10. 5) 

M B 

where energy of midpoint of each energy bin associated with this gain setting 
is E midpoint = M * row number + B (in MeV /3's), 

Columns 21 to 25 (15) Columns 26 to 30 (15) 

Number of rows ignored Number of columns ignored 

Card 4+ Column 2 - namelist format, $FESS_IA=(Number of counts, row 1, column 1); 

then neglecting card column 1, all counts on 64 by 64 analyzer for this gain 
setting, by columns, counts separated by commas 

Card 5+ Column 2 - namelist format, $ FESS_IA=Same but background for this gain 
setting 

Repeat cards 3, 4+, 5+ for other gain settings. If more than one run is to be performed, 
start new set with card 1 and same deck format. 
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APPENDIX B 


PROGRAM LISTING 


C PREJUD PROGRAM MAIN DECK 

COMMON IA 

DIMENSION E( 64,4) , I A( 64,64) , NSUMI (64) ,NSUM2 1 64 ) , MI N R ( 64 ) , M I NB< 6 
IA >, NSB 1( 64) ,NSB2( 64) , XI (64) ,Y1 (64) ,ERl (64) , X? (64) , Y2I64 l,FR2< 64), X 
23 ( 6 4 ) , Y3< 641 ,ER 3< 64) , X4 (64) ,Y4 (64) ,ER4 (64) ,YS1 <6<* I ,YS2< 64) ,YS3< 64) 
3, YS4( 64), XINTI 256) , Y! NT (2561 , YMNT(256) ,ERRO(256 ) ,YTEMP< 256 ), X TEMPI 
4256 ), Z ( 256) ,NR ( 64) 

1 READ! 5,200) 

WR I TE ( 6, 205) 

WRITE! 6, 200) 

READ! 5,202 )N0FD 
IF< N0FD.GT.4) GO TO 1000 
GO TO 100 5 
1000 WP l TE( fc, 201) 

STOP 

ICO 5 DO 3° 0 C I = l.NOFD 

READ! 5,215)XM,B,ICHN, 1 ZE 
215 FORMAT ( 2F 1 0. 5 , I 5 , I 5 ) 

DO 101C J = 1 » ICHN 
1010 E ( J , I ) = 0.0 

ISTAR1 = ICHN + 1 
DO 1015 J = I STAR 1,64 
1015 E< J, I ) = B + XM*F LOAT ( J) 

CALL RD 
M I = 0 

DO 1025 I 1 = I STAR1 ,64 
DO 10 20 J 1 = 1,64 
1020 NR ( J 1 ) = IA( I 1 , J1 ) 

CALL D SUM ( NR , N SUM 1 , NS UM2 , M I , I 1 , 1 ZE ) 

1025 M INRI ID = MI 

C AT THIS POINT IN COMPUTATIONS THE MINIMUM COLUMN NOS. ARE STORED IN 
C. M INRI 64) FOR THIS GAIN SETTING. 

CALL RD 
MR = 0 

DO 104 5 II = I STAR 1 ,64 
DO 10 4 C J 1 = 1 , 64 
1 C4C NR I J 1 ) = I A I I 1 , J 1 ) 

CAL L D SUM (NR, N SB 1, NSB 2, MR, II , 1 Z E ) 

104 5 M IN B( I 1 ) = MR 

C. AT THIS POINT IN COMPUTATIONS THE MINIMUM COLUMN NOS. OF THE BACK- 

C GROUND ARE STORED IN MINS (64) FOR THIS GAIN SETTING. 

GO TO (2000,2200,2400,2600) , I 
C F IRST SET OF DATA 
200 0 DO 200 5 1 2= 1,64 
X II 12 ) = E I 12, I ) 

Y 1 ( 12) = N SUM 2(12 ) -NSB 2(12) 

20C 5 ER1( 121 = SQR T(ABS(Y1 (12) ) ) 

DO 20IC I 3 = I STAR 1 ,64 
IZ = I 3 - ICHN 
X 1 ( IZ > = XI ( I 3 ) 

Yl( IZ ) = Vl< I 3) 



2010 ER 1( i Z ) = ER 1 ( 1 31 

201 1 ICLEAN = 65 - ICHN 

00 2012 I CL E = I C LE AN * 64 

X 1 ( ICLE) =0.0 

Y 1 < ICLE ) = 0.0 

2012 ERH ICLE) = 0.0 

4 FORMAT! 1H1,6X,4 2HFIRST SET OF DATA E NE RGY , COUNTS , AND ERROR.) 

WR ITE! 6,4) 

WR I TE ( 6, 3) (XI ( IK) , Y1 ( I K) ,ER1 ( IK) ,IK = 1 ,64) 

3 FORMAT! 3G 16.8) 

2015 IF( Y 1 ( D.GT.O.O) GO TO 2030 
DO 2020 I 4 = 1,63 
X 1 ( 14) = X 1 ( 14+1) 

Y 1 ( 14) = Y 1 ( 14+1) 

20 20 FR 1( 14 ) = ER.1 ( 14 + 1 ( 

GO TO 201 5 
2C30 NO F 1 = 0 

DO 2035 15 = 1,64 

IF ( Y 1 ( 15) .GT.O.O) N0F1 = NOF 1 + 1 
IF( Yl( 15) .LE.O.O) GO TO 2040 
2035 CONTINUE 
2040 NP = NOF 1 + 1 

DO 2045 I 6 = NP ,64 
X 1( 16 ) = 0.0 
Y 1 ( 16) = 0.0 
2045 FR 1 ( I 6 ) = 0.0 
GO TO 200 0 

2200 DO 2205 I 12 = 1,64 
X 2( 1 1 2 > = E( I 12,1 ) 

Y 2< 112) = NSUM 2(112) - NSB2U12) 

2205 FR 2( 112) = S3R T ( AB S < Y2 ( l l 2 ) ) > 

DO 2210 I 13=1 STAR1 ,64 
IZ 1 = 113 -ICHN 
X 2( IZ 1 ) = X2( I 13) 

Y2< IZ 1 ) = Y2< I 13) 

2210 FR2( IZ 1 ) = ER2( I 13) 

221 1 ICL EAN = 65 - ICHN 

DO 2212 ICLE = ICLEAN,64 
X 2( ICLE) = 0.0 

Y 2 ( ICLE) = 0.0 
2212 ER 2( ICLE ) =0.0 

5 FORMAT (1H1, 6X, 4 3H SECOND SET OF DATA E NERG Y , C OU NT S , AND ERROR.) 
WRITE! 6,5 ) 

WRITE! 6,3 )(X2< IK) , Y2 (I K) ,ER2< I K) ,IK = 1,64) 

2215 IF! Y2< D.GT.O.O) GO TO 2230 
DO 2220 114 = 1,63 
X 2! 114) = X2( 1 14 + 1) 

Y2( I 14 > = Y 2 ! I 1 4 + 1) 

2220 ER2! I 14) = ER2! I 14 + 1 ) 

GO TO 221 5 
2230 NOF 2 = 0 

DO 2235 115 = 1,64 

IF! Y2( 115) .GT.O.O) N0F2 = N0F2 + 1 
IF! Y 2 ( 1155 .LE .0.0) GO TO 2240 
2 235 CONTINUE 
2240 NP ! = NOF 2 + 1 

DO 2245 116 = NP1 ,64 
X 2! I 16 ) = 0.0 
Y21 116 ) = 0.0 
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2245 PR 2 C 1161 = 0.0 
GO TO 300 0 

2400 DO 2405 122 = 1,64 
X 3 ( 12 2) = E ( 122,1 i 

Y 3( 122 ) = N S UM ? ( I 2 2 ) - NSB2 ( I 22 > 

2405 ER 3 ( 122) = SORT (ABS<Y3 <1221 I > 

DO 24 1C 123= I STAR 1 j 64 
I Z 2 = 123 - ICHN 
X 3( I Z 2 ) = X3 ! I 2 3 ) 

Y3( IZ 2 > = Y3( 123) 

2410 FR 3! I Z 2 ) = ER3(I23) 

241 1 l CL EAN = 65 - ICHN 

DO 2412 I CLE = ICLEAN.64 
X 3( ICLE) = 0.0 
Y3< ICLE) =0.0 
2412 ER3< ICLE) = 0.0 

6 FORMAT! 1H 1, 6X ,42HTHIRD SET OF DATA ENERGY, COUNTS, AND ERROR.) 
WR ITE< 6,6 ) 

WR ITE( 6,3 )(X3( IK) , Y3(I K) ,ER3( IK) ,IK = 1 ,64) 

2415 IF! Y3( D.GT.O.O) GO TO 2430 
DO 24 2 C 124 = 1,63 
X 3! 124) = X3( 124 + 1) 

Y3( 124 ) = Y 3! I 24 ♦ 1 ) 

2420 ER 3! 124) = ER3U24 + 1) 

GO TO 2415 
2430 NO F 3 = 0 

DO 24 3 5 1 25 = 1,64 
IF! Y3( 125) .GT.O.O) N0F3 = N0F3 + 1 
IF! Y3( 125 ) .LE .0.0) GO TO 2440 
2435 CONTINUE 
2440 NP 2 = NOF 3 + 1 

DO 2445 I 26 = NP2 ,64 
X 3! 126) =0.0 

Y 3! 126) =0.0 
2445 ER 3! 126) = 0.0 

GO TO 2000 

260C DO 2605 132 = 1,64 
X 4! 13 2) = E ( 132,1 ) 

Y 4! 13 2) = N SUM 2 ( I 32 ) - NSB2U32) 

2605 ER4(I32) = SQR T ( AB S ( Y4 ( I 3 2 ) > ) 

DO 26 1C I 33 = I STAR 1 , 64 
I Z 3 = 133 - ICHN 
X 4! IZ 3 ) = X4( I 33) 

Y4< IZ 3 ) = Y4( I 33) 

2610 FR 4 ( I Z 2 > = ER4! I 33) 

261 1 I CL EAN = 65 - ICHN 

00 2612 ICLE = I C LE AN , 64 
X 4! ICLE) = 0.0 
Y 4 ( ICLE ) =0.0 
2612 ER 4! ICLE ) =0.0 

7 FORMAT! 1H 1, 6X , 43HF0UR TH SET OF DATA E NE RG Y , COUNTS , AND ERROR.) 
WR ITE! 6, 7 ) 

WR ITE ( 6, 3 ) ( X4< IK) , Y4! I K) ,ER4 ( I K) ,1 K=i ,64) 

2615 IF! Y4( D.GT.O.O) GO TO 2630 
DO 262C 134= 1,63 
X 4 ( 134) = X4( I 34 + 1) 

¥4(134) = Y4( I 34 «- 1 ) 

2620 ER 4! I 34 ) = F R4! I 34 «■ 1) 

GO TO 2 615 



2 63 C MO F 4 = 0 

DO 2635 I 35 = 1 ,64 
I Ft Y4I 135 ) .GT.O.O) N0F4 = N0F4 + 1 
IF( Y4( 135 I .LE .0.0) GO TO 2640 
2635 CONTINUE 
264 C MP3 = NOF 4 4- 1 

00 2645 I 36 = NP3 ,64 
X 4( 136 ) = 0.0 
V4I 136 > = 0.0 
2645 ER4I 136) =0.0 
GO TO 300 0 
3000 CONTINUE 

N ST = NOF l 

C FERE THE FIRST SET OF DATA IS SMOOTHED 
3010 CALL SMI Y l » NOF 1 » YS1 ) 

C CALL SSMI YSl.NOFl) CARDS MAY BE ADDED OR TAKEN OUT TO FURTHER 
C SMOOTH DATA, OR INCREASE FIDELITY TO ORIGINAL DATA. HERE 2 HAVE 
C BEEN USED. 

CALL SSMI YSl.NOFl) 

CALL SSMI YSl.NOFl) 

IF! NST .EQ .NOF l ) GO TO 3020 
NO = NOF 1 «■ 1 
DO 3015 IM = ND.64 
X II IM ) = 0.0 
Y S 1 1 IM ) = 0.0 
3015 ER1I IM ) = 0.0 

3020 IFUNOFD - D.LE.O) 30 TO 3200 
NST = NOF 2 

C HERE THE SECOND SET OF DATA IS SMOOTHED 
3030 CALL SMI Y2.N0F2.YS2) 

C CALL SSMI YS2.NOF2) CARDS MAY BE ADDED OR TAKEN OUT TO FURTHER 
C SMOOTH DATA, OR INCREASE FIDELITY TO ORIGINAL DATA. HERE 2 HAVE 
C BEEN USED. 

CALL SSMI YS2.NOF2) 

CALL SSMI YS2.NOF2) 

IF! NST.EQ .N0F2) GO TO 3040 
ND = NOF 2 + 1 
DO 3035 IM = ND , 64 
X2I IM ) = 0.0 
YS2I IM I = 0.0 
3035 ER2I IM ) = 0.0 

3Q4C IFUNOFD - 2I.LE.0) SO TO 3200 
NST = NOF 3 

C HERE THE THIRD SET OF DATA IS SMOOTHED 
3050 CALL SMI Y 3, NOF 3, YS3) 

C CALL SSMI YS3.NOF3) CAROS MAY BE ADDED OR TAKEN OUT TO FURTHER 
C SMOOTH DATA, OR INCREASE FIDELITY TO ORIGINAL DATA. HERE ? H AV F 
C BEEN USED. 

CALL SSMI YS3.N0F3) 

CALL SSMI YS3.N0F3) 

IF! NST.E3 .N0F3 ) GO TO 3060 
ND = NOF 3 + 1 
DO 3055 IM = ND , 64 
X3I IM ) = 0.0 
YS3I IM ) = 0.0 
3C55 FR 3 1 IM > = 0.0 

3060 I F { INOFD - 3) .IE. OS 30 TO 3200 
NST = NOF 4 
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C HERE THE FOURTH SET OF DATA IS SMOOTHED 
3070 CALL SM ( Y 4 ,NQ F 4 * Y S4 ) 

C CALL S SMI Y S4 , N OF 4) CARDS MAY BP ADDED OR TAKEN OUT TO FURTHER 
C SMOOTH DATA ,0R INCREASE FIDELITY TO ORIGINAL DATA. HERE ? HAVE 
C BEEN USED. 

CALL SSM( YS4.N0F4) 

CALL S$M( Y S4 » NOF 4 ) 

IFINST.EG .N0F4) GO TO 3200 
MD = N0F4 «■ 1 
DO 3075 IM = ND,64 
X 4! IM ) = 0.0 
YS4! IM ) = 0.0 
3075 ER4( IM ) = 0.0 
GO TO 220 0 

3200 GO TO (3300,3400,3600,3800), NOFD 
3300 DO 3305 IPIC = l.NOFI 
X IN T ( IPIC ) = X 1 ( IPIC) 

Y IN T( IPIC ) = YS1 ( IPIC ) 

3205 E R R 0 ( IPIC ) = ER1( IPIC ) 

NALL = NO F 1 
GO TO 400 0 

3400 CSUB2I = CONI XI,YS1 ,N0F1,X2 ,YS2 , NOF 2 ) 

8 FORMA T( IH5, 3X, 53HNORMALI ZA TI ON CONSTANT OF SECOND DATA SET TO FIRS 
XT IS) 

WR ITE( 6,8 ) 

WR I TE ( 6,6010) CSUB21 
DO 34 0 5 IF 21 = 1 , NOE 2 

Y S 2 ( l F 2 1 ) = CSUB21*YS2(IF21 ) 

3405 ER 2 ( I F 2 1 ) = C SUB 2 1*FR 2 ( I F 2 1 ) 

0 FORMA T ( 5X , 26HN0RMAI. I ZED SECOND DATA SET) 

WRITE! 6,9) 

WRITE! 6,6010) ( X2 ( I ) , YS2 ( I ) ,E R2 ( I ) , I = 1,64) 

DO 3410 IRIC = 1.N0F1 
X IN T ( IRIC ) = XI (IRIC) 

YINT! IRIC ) = YSl(IRIC) 

34 1 C ERRO! IRIC ) = E R 1 ( IRIC ) 

CALL MERGE! XINT, YINT, ERRO, NOF 1 , X2 , YS2 , ER2 , NOE 2 , NALL ) 

10 FORMAT! 1H1,5X,33HMER3FD FIRST AND SECOND DATA SETS) 

WR ITE! 6, 10) 

WR I TE! 6, 6010) ( XI NT ( I ) , YI N T ( I ) , E RRO ( I ) , I = 1,130) 

GO TO 400 0 

360 C CSUB21 = CON! XI ,YSl ,N0F1 , X2 ,YS2 ,N0F2) 

WRITE! 6,8 ) 

WRITE! 6,6010) CSUB21 

CSUB32 = CON! X2,YS2,NOF2,X3,YS3, NOF 3 ) 

11 FORMA T ( 4X , 53HNORMAL I Z ATI ON CONSTANT OF THIRD DATA SET TO SECOND IS 
X ) 

WR ITE! 6, 1 1 ) 

WRITE! 6,6010) C SUB 32 
DO 3605 I F 2 1 = 1.N0F2 
YS2UF21) = CSUB21*YS2(IF21 ) 

360 5 E P 2 ( I F 2 1 ) = C $UB21*ER2 ( I F 2 1 ) 

WR ITE! 6,9 ) 

WRI TE! 6,6010) ( X2( I ) ,YS2( I ) ,ER2( I ) ,1=1 ,64) 

DO 36 10 IF 31 = 1 , NOF 3 

Y S 3! I F 2 1 ) = C SUB32*CSUB21*YS3 (IF31 ) 

3610 ER 3( I F 2 1 ) = C $ UB3 2*C SUB 2 1 *E R 3 ( I F 31 ) 

1? FORMAT! 5X , 25HNOPMAL I ZED THIRD DATA SET) 

WR I TF ( 6, 1 2 ! 
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WR I T e ( 6,6010 ) { X 3 ( I ) ,YS3(I ) ,ER3 ( I > * I = 1 ,64) 

DO 36 15 I SIC = 1 ,N0F1 
X IN T( I SIC ) = XI ( I SIC I 

Y INTI I SIC ) = YS1( I SIC) 

3615 ERRO! I SIC ) = ERK I S I C i 

CALL MERGE! X I N T , Y I NT , FRRO , NOF 1 ,X2,YS2 ,ER2 ,N0F2 , NALL) 

WR ITEI 6, 10) 

WR I TE( 6,6010) ( XI NT (I ) ,YlNT(l ) , ERROI I ) ,1=1 ,128) 

N0FM1 = NALL 

CALL MERGE ( XINT, Yi NT,ERRO, NOF Ml , X3 , YS3 ,ER3 ,N0F3 , NALL ) 

13 FORMAT! 1H 5 , 6X , 3 8HN0RM A L I Z A TI ON OF FIRST THREE DATA SETS) 

WR ITE! 6, 1 3) 

WP. I TE( 6,6010) (XI NT (I ) ,YINT(I ) , ERR 0(1 ) ,1 = 1 ,196) 

36 20 GO TO 400 0 

3800 C SUB 2 1 = CON! X 1 , YS1 »N0F1 ,X2 »YS2* NOF 2 ) 

WR I TE( 6,8 ) 

WRITE! 6,6010) C SUB 21 

CSUB32 = CON! X2,YS2 ,NOF2,X3 ,YS3 ,N0F3) 

WRI TE( 6, 1 1 ) 

WRITE! 6,6010) CSUB32 

CSUB43 = CON! X3,YS3»N0F3 , X4,YS4 , NOF 4 ) 

14 FDRMAT(4X» 53HN0RMAL I ZA TI ON CONSTANT OF FOURTH DATA SET TO THIRD IS 
X) 

WR ITE! 6, 14) 

WRI TE< 6,6010) C SUB43 
DO 3805 I F 2 1 = 1.N0F2 

Y S 2! I F 2 1 ) = CSUB21 *YS2!IF21) 

3805 ER 2! I F 2 1 ) = CSUB21 * ER2UF21) 

WR ITE! 6,9 ) 

WR ITE! 6,6010 ) ( X2< I ) ,YS2( I ) ,ER2( I ) , I = 1,64) 

DO 38 1 C I F 31 = 1 , NOF 3 

YS3IIF31) = CSLB32*CSUB21*YS3 ( I F 31 ) 

3810 ER 3! IF 31 ) = C SUB3 2*C SUB2 1 *E R 3 ( I F 31 ) 

WR ITE! 6, 12) 

WRITE! 6.6010HX3! I) ,YS3(I ) ,ER3(I ),I = 1,64) 

DO 3815 I F 41 = 1 , NOF 4 

YS4! IF 41) = CSUB43*CSUB32*CSUB21*YS4(IF41 ) 

3815 ER 4 ( I F 4 1 ) = C SUB43*CSUB32*C SUB2 1 *ER4 ( I F41 ) 

15 FORMAT ( 5X , 26HN0RMALI ZED FOURTH DATA SET) 

WP ITE! 6, 15) 

WRITE! 6»6C10) (X4( I) ,YS4(I ) ,ER4(I ) , I = 1,64) 

DO 3820 I TIC = 1 ,N0F1 
X IN T< I TIC ) = XI ( I TIC ) 

Y INT! I TIC ) = YS1 ! I TIC ) 

3820 ERROI ITIC ) = ER 1 ! I TIC ) 

CALL MERGE ( X I N T , Y I NT , ERRO , NOF 1 , X2 , YS2 , ER2 , NOF 2 , NALL ) 

WR ITE ( 6, 10) 

WR I TE( 6,6010 ) ( XINT ( I ) , YINT! I ) ,ERRO( I ) , I = 1,128) 

N0FM1 = NALL 

CALL MFRGE! XINT, YINT, ERRO, NOE Ml , X3 , YS3 * E R? , N0F3 , NALL ) 

WR ITE! 6, 13) 

WR. I TE! 6,6010) ( XINT! I ) , YINT! I ) ,E RRO! I ) , I = 1,192) 

N0FM2 = NALL 

CALL MERGE ( XINT , YINT, ERRO, NOF M2 , X4 ,YS4 ,ER4 ,N0F4 , NAl L) 

16 FORMAT! 1H5,6X, 3 7HN0RMA L l ZATI ON OF FIRST FOUR DATA SETS) 

WR ITE! 6, 16) 

WR I TE ( 6,6 010) ( X INT! I ) , YI NT! I ) ,E RRO ( I ) , I = 1 , 256) 

3825 GO TO 4000 
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C EXPONENTIAL SMOOTHING OF DATA POINTS 

4000 NALLM = NALL - I 

00 40 1C f ZR = 2, NALLM 


D1MEN S ION 

E XSM ( 3) ,E YS M ( 3 ) ,EWSM(3) 

YMN T( l ) 


Y I N T ( 1) 



EX $M{ 1 ) 

= 

XINKIZR - 

1 ! 


6XSM( 2 ) 


X IN T ( IZR! 



EX SMI 3 ) 

= 

XINT ( I ZR 4- 

1) 


F.Y SM { 1 ) 

= 

YINTI I ZR - 

1 ) 


EY SM( 2 > 

= 

Y INK IZR) 



EY SM( 3 ) 


YINTI I ZR <• 

1 ) 


FWSMI 1 ) 

= 

FRR 0 ( I ZR) /ERRO! I ZR 

- 1) 

EWSM( 2 ) 


1.0 



EWSM( 3 ) 

= 

ERROd ZR) /ERROI I ZR 

4- 1 ) 

4C10 YMNT! I ZR ( 

= WEXSMO (EXSM,EYSM, 

EWSM) 


YMN T(NALL ) = YINT (NALL) 

NALL = NALL + 1 

C INTERPOLATION POINTS FOR FEROOR 

DAT A { l ( I ) , 1 = 1 ,256) /. 0317,. 0381 ,.0444, .0508 , .0572 , . 06? 4, .3762, .0-488 
X, . 1015,. 1 20 5,. 1395, .165,. 196 8, .22 86 , .2 64 7,. 3048,. 3 49 2, . 39 37, .4 44 5, 

X . 49 53 , .55 24, . 6096, .6668,. 7239,. 7810 ,.8332 ,.895 4,. 9 525, .95 25, 1.016, 

XI. 0 16, 1.080, 1.080, 1.143, 1.143, 1.206, 1.27, 1.334, 1.397,1. 46, 1.52 4, 1 . 
X5875, 1 .65 1, 1. 714, 1. 778,1. 842 , 1 . 905 , 1 . 968 ,2.9 32 ,2.096, 2 . 159. 2. 222, 2 
X. 286, 2. 35 0,2. 413, 2. 476, 2. 54, 2. 66 7, 2. 794, 2.92 L ,3.043,3.175, 3.302, 3. 
X4 29, 3.556, 3. 6 8 3,3. 81,3.937,4. 064 ,4.191 ,4 . 3 1 8 , 4 . 445 , 4 . 5 7 2 , 4 . 699 , 4 . 3 
X26, 4. 953, 5. 08, 5. 20 7, 5. 334, 5. 44, 5. 58, 5. 72, 5. 84, 5. 96, 6. 09, 6. 2 2, 6. 34, 
X6. 5 4, 6. 73, 6. 92, 7. 11, 7. 30, 7. 49, 7. 68 ,7. 88, 8. 07, 8. 26, 8. 44, 8. 6 4, 8. 82, 9 

X. 02, 9. 21, 9. 40, 9.58, 9. 78,9. 97, 10. 16, 10. 35 , 10. 55, 10. 73, 10. 93, 11. 10,1 

XI. 30, 14340.0/ 

C SPLINE FIT 

5000 NALL = NALL - l 

CALL SPLN ( XI NT, YMNT ,NALL» Z.YTEMP ,1 13) 

CALL ERROR ( XI NT, ERR 0, NALL , Z , X TE Mp, 1 13 > 

5010 KAS = 1 

5015 [F( YTEMP ( KAS) .NE.0.05 GO TO 5030 
C LOW ENERGY EXTRAPOLATION 

5020 CALL EXTRA! XINT(i) , YMNTU ) ,XI NT ( 2) , YMNT ( 2 ) , X I NT ( 3 ) , Y MNT ( 3 ) , Z ( < AS ) , 
1Y TEMP ( KAS ) ,KA S) 

XTEMP ( KAS ! = Y TEMPI KAS)* (ERROd ) / YMNT ( 1 ) ) 

KAS = KAS «• 1 
GO TO 5015 
5030 N IR = 1 

DO 5035 L INC = 1,256 
I F( YMN T(L INC ) .GT.25.0) NALL = L I NC 
IF< YMNT(L INC ) „LE. 25.0) GO TO 5036 

5035 CONTINUE 

5036 CONTINUE 

DO 504C LR = 1,113 
IFIYTEMP(LR). GT.25.0) GO TO 5040 
C HIGH ENERGY EXTRAPOLATION 

CALL HIEXTI XINKNALL-4) , YMN T ( NA L L-4 ) , XI NT (NALL-3) , Y MNT ( N AL L - 3 ) , 

IX IN T( NALL- 2) , YMNT(NALL-2) , XI NT (NALL-1 ) , YMNT (NALL-1 I ,X I NT ( N ALL ) > Y M N 
2T(\ALL ),Z ( LR ) , YTEMP (LR) ,NIR) 

NIR = 2 

XTEMP ( LR ) = ERRO (NALL) 

504C CONTINUE 

6 COO WRITE! 6 , 6 010 S ( Z ( 1 ) , YTE MP ( I ) ,XTE MPH ) ,1=1 ,113) 

200 FORMAT ( BOH 

1 ) 
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20] FORMA T ( 19 H TOO MANY DATA SETS) 

202 FORMAT] 15 ) 

205 FORMAT] 1H 1 ) 

6 C 1 0 FORMAT] 30 16.8) 

6^20 PUNCH 200 

PUNCH C03 0,( I , YTEMP] I ) .XTEMPU) ,1 = 1 ,1 13) 
6030 FORMAT] 25X , I 3 , 22X ,E 1 O'. 3 , 1 OX ,E 1 0. 3) 

GO TO 1 
FND 


$ I EFTC L FNER DECK 

SUBROUTINE ERR OR ( X, Y , N , Z , XT , I NT I 
C ERROR ESTIMATE 

DIMENSION X< 256 ),Y( 256) , Z ( 2 56 ) , XT! 2 56) 

DO 40 I = 1, INT 

IF] Z ( I ).L T.XI 1 ) .OR . Z( I ) .GT. XI N) ) XT ( I ) = 0.0 

IF] Z< I ).L T.X( 1 ) .OR. Z(I).GT. X< N) ) GO TO 40 

IF] Z ( I ).EQ.X< I) ) XT { I ) = Y ( 1 ) 

IF] Z( I ).EQ .X( 1 ) ) GO TO 40 
I F ( Z ( I ).EQ.X(N> ) XT ( I ) = Y ( N) 

I F( Z( I ).EQ .X(N ) ) GO TO 40 
DO 30 J = 2, N 

I F ( Z ( I ) .1 E .X ( J ) . AND . Z ( I ) . GE . X ( J-l ) ) XT(I) = Y(J-I) +( ( Y ( J ) -Y ( J - 1 ) ) 
1/ ( X(J ) — X ( J — 1 ) ) )*( Z( I)-X(J-i) ) 

30 CONTINUE 
40 CONTINUE 
RETURN 
END 


SieFTC WEXS DECK 

FUNCTION WEX SMO ( X , Y , W ) 

C WEIGHTED EXPONENTIAL SMOOTHING OF DATA 
DIMENSION X(3),Y(3) ,W(3) ,H(3) 

DOUBLE PRECISION ZUCK , ZAC K , ZORP , AKE , AKK , AR P, F ARGO, A PPL E 
NAG = C 
DO II K = 1,3 
IF] Y( K ) >1 2,11 ,11 
12 Y(K ) = A3 S( Y ( K ) ) 

NAG = 1 

11 CONTINUE 

I F( NAG ) 13,60,13 

12 WRITE ( 6,6 ) 

6 F0RMAT(77HJNEGATI VE Y VALUES HAVE BEEN MADE POSITIVE IN CALCULATI3 
INS BEYOND THIS LINE. Ill S 
60 DO 130 I =1,3 
H( I ) = AL 0 G { Y ( I ) ) 

130 CONTINUE 

ZUCK = 0. ODO 
ZACK = 0.000 
ZORP = 0.000 
AKE = O.ODO 
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AKK = O.ODO 

ARP = 0.000 

DO 10 I = 1,3 

ZUCK = ZUCK + Will 

ZACK = ZACK + X( I )*W( I ) 

zorp = zorp + h ( n + wm 

AKE = AKE t K(I )*W( I ) 

AKK = AKK + W( I )* X( I )**2 

10 ARP = ARP 4- XI I )*H ( I )*W(I I 
FARGO = Z UCK* AKK - ZACK*AKE 
IF( FAR GO ) 30, 20,30 
20 WRITFI 6,8) 

8 FORMA T ( 30H DIVISION BY ZERO NO SOLUTION.) 
RETURN 

30 APPLE = ( ZORP*AKK - ZACK* ARP) /FARGO 
ALPHA = DEXP ( APPLE ) 

BETA = ( Z UCK* ARP - ZORP* AKE ) /F ARGO 
WEXSMO = ALPHA*E XP ( BE TA* X (2 ) ) 

RETURN 

END 


$ I BFTC NOW DECK 

SUBROUTINE SSM(X,N) 

C FIVE POINT SMOOTHING OF DATA POINTS 

DIMENSION X( 64 ) 

DOUBLE PRECISION XD ( 64 ) , X L ( 64 ) 

DO 10 I = 1, N 
10 XD( I) = XI I) 

IFIN.LT.5) GO TO 50 
DO 20 I = 1 , N 

IF< XD( D.GT.O.O) XL ( I ) = D L0G10 ( XD ( I ) ) 

20 CONTINUF 

X D ( 1 ) = ( XL ( 1 )*69. ODO 4- XL(2)*4.0D0 - Xl(3)*6.0DP 4- XL(4l*4.0DO - 

1XL( 5) ) /70.0D0 

XD ( 2 ) = ( XL ( 1 )* 2. ODO 4- XL(2)*27.0D0 4- XL ( 3 ) * 1 2 . 0 DO - X L ( 4 ) *8 . 0 DO 4- 

1 XL ( 5 )*2. ODO) /35. ODO 

XD(N-l) = ( XL ( N-4 > * 2. ODO - XL I N-3 ) *8 . 0 DO * X L ( N-2 ) * 1 2 . 0 DO + XLIN-1 
11*27. ODO + XL (N )*2. ODO) /35. ODO 

XO(N ) ={-XL (N-4) + XLI N-3) *4. ODO - XL(N-2)*6.0D0 4- X L< N- 1 ) * 4 .HDD 4- 
1XL(N)*69.0D0) /70.OD0 
NM 2 = N - 2 
DO 30 I = 3, NM2 

30 XD( I ) = ( — XL ( I-2)*3.0DO f XL I I - 1 ) * 1 2 . ODO 4- XL(I)*17.0DO 4- XL(l4-l>* 

112. ODO - XLI I 4- 2 ) * 3 . 03 0 ) /3 5 . OD 0 
DO 40 I = 1 , N 
40 X( II = I 1 C.0D0)**XD (I ) 

5C RETURN 
END 
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sief fc nun 


DECK 


SUBROUTINE SM ( X , N * X S) 

C FIVE POINT SMOOTHING OF OATA POINTS 

DIMENSION! X( 64) ,XS( 64) 

DOUBLE PRECISION XD ( 64 ) , X L ( 64 S 
DO 10 I = 1 1 N 
10 XD( I) = X ( I ) 

IF(N.LT.S) GO TO 50 
DO 20 f = I * N 

IF( XD( I I.GT.O.O) XL ( I I = D LOG 1 0 ( XD ( I ) > 

20 CONTINUE 

X D ( 1 ) = ( XL( I (*69.000 4 XL<2 I $4. ODO - XL(3)*6.OD0 + XL ( 4 ) *4 , n D0 - 

IXL ( 5 ) ) /70.000 

XD ( 2 ) = I XL ( 1 ) * 2. 00 0 4- XL(2)*27.0D0 4- XL(3>*12.000 - XL(4)*8.0D0 4 
I XL ( 5 )*2. ODO ) / 3 5. 00 0 

XD(N-l) = ( XL (N-4)*2. ODO - XL ( N-3 ) *8 . 0 DO 4 X L ( N-2 ) * L 2 . 0 DO 4 XL C N-l 

11*27. ODO 4 XL(N)*2.0D0) /35.000 

XO(N ) = ( - XL (N- 4 ) 4 XL ( N-3 ) *4. ODO - XL I N-2 > *6 . 0 DO 4 X L ( N- I ) *4 .ODO 4 
IXL ( N ) * 69. 000) / 70. ODO 
NM 2 = N - 2 

DO 30 I = 3» NM2 

30 XD( I) = ( -XL ( 1-21*3.000 4 XL ( I - 1 ) *1 2 . 0 DO 4 XL(I)*I7.0Dn 4 XL(l4l)* 
112. ODO - XL( I42)*3.000) /35.0D0 
DO 40 I = 1 1 N 
40 XS( I) = I 10. ODO )**XD( I ) 

R E T UP N 

50 DO 60 I = I ? N 
6C X S ( I) = X ( I ) 

RETURN 

END 


SieFTC MUB DECK 

SUBROUT IN E MERGE ( XINT S VINT ,ERRO, NOF 1 » X2 ? Y2 ? E R2 ? N0F2, IF INAL ) 

C MERGE ROUTINE 

D IMEN S ION XI ( 2 56) »Y1 ( 2 56) S ER1 (2 56) , X2 (64) , Y2 < 64 I , E R2( 6 4 ) , X INT ( 256 ) 
I » Y INT ( 256 ) fERRO ( 2 56 ) 

IF(N0F1.GT. 256) STOP 
DO 10 I = 1 » NOF 1 
X 1 ( I ) = XINTt I ) 
v 1 ( I ) = V INT( I ) 

ERKII = ERRO ( I ) 

X IN T( I ) = 0.0 
Y IN T( I ) = 0.0 
1 C F R P 0 ( I ) = 0.0 
15 IFF = 1 
IMF = 1 

[FINAL = 1 

20 I F < X 1 ( IFF >.LT.X2( IMF) ) GO TO 45 
I F ( X 1 ( IFF ).F0.X2< IMF) ) GO TO 65 
25 X INK I F IN AL ) = X2 (I MF ) 

YINK I FINAL! = Y2 ( I MF ) 

ERR 0 ( IF IN AL ) = ER2< IMF ) 

30 !F( IMF.EO .N0F2) GO TO 35 
IMF = IMF 4 I 
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{FINAL = IF INAL + 1 
GO TO 20 

35 I F ( IFF.E3.N0F1) GO TO 40 
IFINAL = {FINAL + I 

X IN T( {FINAL) = XI (IFF ) 

Y IN T( IFINAL) = Yl ( IFF ) 
ERROI I F IN AL > = FR 1 ( IFF ) 

IFF = IFF 4- I 
GO TO 35 


4 0 


44 

45 


50 


55 


60 


65 


IFINAL = IFINAL 
X INK I FINAL ) = 

Y IN T( IFINAL ) = 
ERRO< IFINAL) = 
RETURN 

X INTI IFINAL ) = 
YINTI IFINAL) = 
ERROI IF IN AL ) = 

I F ( IFF .F3 .N3F I) 
IFF = IFF + 1 


+ I 

XI (NOF 1 ) 
Yl I NOF 1 ) 
ERKN0F1 ) 

XI ( IFF ) 

Yl ( IFF ) 
FRMIFFt 
GO TO 55 


IFINAL = IFINAL + 1 

GO TO 20 

I FI IMF.E3 oNO F 2 ) GO TO 60 
IFINAL = IFINAL + I 
X IN T I IFINAL) = X2IIMF ) 

YINTI IFINAL) = Y2 1 1 MF ) 

ERROI IFINAL ) = ER2I INF ) 

IMF = IMF 4- l 
GO TO 55 

IFINAL = IFINAL +■ 1 
X INTI IF IN AL ) = X2 (NOF 2 ) 

YINTI IFINAL) = V2INOF2) 

ERROI I FINAL) = ER2IN0F2) 

RETURN 

IF( ER II IFF ) .LT.ER2I IMF ) ) GO TO 70 
XINTI IFINAL ) = X 2 1 1 M c ) 

Y INTI IFINAL ) = Y2 I I MF ) 

ERROI I FIN AL ) = ER2 I I MF ) 

IF( IFF .F3 .NOF l ) GO TO 44 
IFF = IFF 4- 1 


GO TO 20 

70 X INTI I FINAL ) = XI I IFF ) 

YINT( IFINAL) = Yl (IFF ) 

ERROI IFINAL) = FRl(IFF) 

IF! IMF . f 3 .NOF 2 ) IFF = IFF 4- 1 
IF( IMF.E3 .NOF 2 ) GO TO 35 
IMF = IMF 4- 1 
GO TO 50 
END 


$ IBFTC PR EX TR DECK 

SUBROUTINE E X TR A I X 1 , Yl , X2 , Y2 , X3 , Y3 , XI NT , Y I NT , KR2 > 

C LOW ENERGY EXTRAPOLATION 

DIMENSION XI 3 ) , YI 3) ,WI 3) ,H I 3) 

DOUBLE PRECISION ZUCK , ZAC K , ZORP , AKt , AKK , ARP s F ARGO, A PP L F 
IF ( KR 2 .GT . 1 >G0 TO 100 0 



NAG = C 

W( 1 ) = 10.0 

W ( 2 ) = 1.0 

W( 3 ) =0.10 

X ( 1 ) = XI 

X< 2 ) = X2 

X ( 3 ) = X3 

Y( 1) = Y 1 

Y( 2 ) = Y2 

Y { 3 ) = Y3 

DO 11 K = 1,3 

IF! Y( K ) )1 2, 11 , 11 

12 Y ( K ) = Ad S( Y( K ) ) 

NAG = 1 

1 1 CONTINUE 

IF(MAG) 13,60,13 

13 WRITE! 6,6) 

6 FORMAT! 77 HJN EGA TI VE Y VALUES HAVE BEEN MADE POSITIVE IN CALCILATIO 
INS BEYOND THIS LINE. ///> 

6C on 130 I = 1,3 

H( I ) = ALC1G! Y! I ) ) 

130 CONTINUE 

ZUCK = 0.000 
ZACK = 0.000 
ZORP = 0.000 
AKE = 0.000 
AKK = 0.000 
ARP = 0.000 
DO 10 I = 1,3 
ZUCK = ZUCK + W(I) 

ZACK = ZACK + X! I )*W( I ) 

ZORP = ZORP + H ( I )*W( I ) 

AKE = AKF + X! I >*W(I) 

AKK = AKK + W ( I ) * X ( I ) ##2 

10 ARP = ARP + X! I )*H( I )*W(I ) 

FARGO = Z UCK* AKK - ZAC K*A KE 
I F ( FAR GO ) 30, 20,30 

20 WRITE! 6,8) 

8 FORMA T ( 30H DIVISION BY ZERO NO SOLUTION.) 

R F T UR N 

30 APPLE = ( ZORP*AKK - ZACK# ARP) /FARGO 
AL°HA = DEXP (APPLE) 

BETA = ( Z UCK*ARP - ZOR P*A KE ) /F A RGO 
1000 YINT = ALPHA * E XP ( BE TA* X I NT ) 

RETURN 

END 


SIBFTC UP DECK 

SUBROUTINE DSUM (NROW, NSUM1 ,NSUM2 ,MI * I , I Z E ) 

C MINIMUM FINDING ROUTINE 

DIMEN S ION NROW! 64) , NSUM1 (64) , NSUM2 (64 ) ,A(6A ) 
MUST = MI 

2 FORMAT! 1H , 9HDSUM1 ERR, I 6) 

K = I 

WRITE! 6, 3 INROW, I 
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3 FORMAT (161 7/1 61 7/161 7/1 61 7/1 5) 

Of) 10 II = 1 , 64 
10 A < I I ) = NROWI I I ) 

IF! MI .EQ.OI GO TQ 12 
IF(MI.EO.I?El GO TO 60 
GO TO 100 

12 IF! A! I ZE+ 1 ! .GT. 100. O.OR.Ad ZE ♦ 2 > .GT. 130.0} GO TO 20 
I F( A( IZE + 3) .GT. 100. 0. OR. A ( I ZE«-A) .3 T. 100.0 ) GO TO 20 
GO TO 60 
20 IA = 2 + I ZE 

25 I F( A( IA-1 ) .GT.A! I A ) .AND. A ( I A) .LT. A! I A* 1 ) ) GO TO 3D 

IF! A! I A- 2 ) .GT.A! I A) .AND. ( A ( I A) . E 0. A (I A-l) . AND. A! I A ) . LT . A( I A+ 1 ) ) ) G 
10 TO DC 

IF! ( ( A! IA-3) .GT.A! I A) . AND. A! I A) .EQ. A( I A-l ) ). AND. A! I A) . FQ. A! I A- 2) ) . 
1AND.A! IA) .IT. A ( IA+1)) GO TO 30 
IF! A! I A ) . E Q . 0 . 0 ) GO TO AO 

26 IA = IA-*- 1 

IF! IA.CE. 6A1WRI TE (6,2) IA 
GO TO 25 

30 IF! A! IA) . GE.( A! IA + D-3. C* SQR T ( A ( I A* 1 ) ) ).OR.(A( I A) .GT.A! IA+2) ) ) GO 

XTO 35 
GO TO AO 

35 IF! (A! IA+2)-A( IA) I.LT.C.OIGO TO 26 
IF! (A! IA+3)-A ( IA) ).LT. 0.0)50 TO 26 
IF! (A! IA+AJ-A! IA) ).LT.O.O)GO TO 26 
AO I SUM 1 = 0 
I SUM 2 = 0 
IB = IA-1 
IC = IA *■ 1 
DO A5 IR = 1, IB 
A 5 I SUM l = I SUM1 + NROW! I R) 

I SUM 1 = I SUM 1 «- NROW! I A ) / 2 

I SUM 2 = NROW! I A ) / 2 
DO 50 IS = IC.6A 

50 ISUM2 = I SUM 2 + NROW! IS) 

I F( NRO W( IAl.GT. ISUM2/ 7) I SUM2 = 0 
IF(NR0W(IA).GT.ISUM2/ 7) I SUM1 = 0 
I F ( NR 0 W( I A ) . G T . I SUM 2/ 7) I A = 0 
IF! ISUM2 . EQ . 0) GO TO 51 

IF! ( ISIM1/ISUM21.LE. 10000) GO TO 51 
IA = 0 
I SUM 1 = 0 
I SUM 2 = 0 

5 1 M I = I A 

N SUM 1 < K > = ISUM1 
N SUM 2 ( K ) = ISUM2 

7000 FORMA T ( 5X, 2AHMINI MUM COLUMN ABOVE IS ,12) 

WRITE! 6 , 7 000 ) M I 
RETURN 

60 IF! AM AXi! A! IZE + l) ,A (I ZE + 2) ,A( I ZE+3) ) .EQ.O .0) G0T065 
GO TO 1 5 
65 IA = IZF 
I SUM 1 = 0 
I SUM 2 =0 
DO 70 IZ = 1 * 6 A 
70 I SUM 2 = I SUM 2 + NROW! I Z) 

GO TO 51 
75 IA = IZE 1 3 
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76 IF< A( IA).EQ.O.O.AND.A ( IA-1 >.EQ.O.O)GO TO 80 
IA = [A +• 1 

! F ( I A .EQ. 64) WR I TE f 6,2) I A 
GO TO 76 
8 0 I A = I A- 1 

I SUM 1 = 0 

I SUM 2 =0 

DO 85 IZZ = 1 , I A 
85 I SUM 1 = I SUM 1 NROW( I ZZ) 

03 90 IZZ 1 = I A , 64 
90 I SUM 2 = I SUM 2 + NR0WUZZ1) 

G3 TO 51 

100 IFIMI-5.GT.IZE) IA =MI - 5 
I F ( M I-5.GT.IZE ) GO TO 25 
GO TO 20 
END 


% IBFTC Z ING DECK 

FUNCTION CONI X,Y,N,XF ,YF,M) 

C NORMALIZATION CONSTANT ROUTINE 

D IMEN S ION XI 2 56) , Y< 25 6) ,X1 (256) , Y1 (256) ,YS (256 ) ,XF( 64 ) , YF( 64) 
DO 10 I = 1,64 
X L ( I) = X F ( l ) 

Y 1 ( I > = Y F ( I ) 

1 C Y S ( I ) =0.0 

CALL SPLN (X,Y,N,XL,YS,M) 

CON F =0.0 
DIV = C.O 
DO 20 I = 1,64 

IF( YS( D.NE.O.OIDIV = DIV+ 1.0 

IF< YS< U.NE.0.0) CONF = CONF 4- YS(I)/Y1<!> 

20 CONTINUE 

CON = CONF /DIV 

RETURN 

END 


tIBETC DUB DECK 

SUBROUTINE RD 
C READ ROUTINE 

COMMON IA 

DIMENSION I A ( 64 , 64 ) 
NAME! I ST/FESS/IA 
1 READ! 5, EE SS) 

RETURN 

END 
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$ IBFTC POSTEX DECK 


SUBROUT IN E H I E X T { X 1 » Y l » X2 , Y2 , X3 , Y3 * X4 , Y4 , X5 , Y5 , X I NT , Y I NT , K 8 8 ) 

C HIGH ENERGY EXTRAPOLATION ROUTINE 
DIMENSION X( 5) , Y( 5) »W( 5) ,H< 5) 

OOUBL E P 3 EC I SI ON ZUCK , ZAC K ,ZORP , AKE , AKK. , ARP , F ARGO, A PP L E 

IFIK88.GT .1) GO TO IDOO 


NAG = 


C 

W( 1 ) 

= 

1.0 

W( 2 ) 

= 

2.0 

W( 3 ) 

= 

4.0 

W( 4 ) 

= 

8. 0 

W( 5 ) 

= 

16.00 

X! 1 ) 

= 

XI 

X! 2 > 

= 

X2 

X! 3 ) 

5S 

X3 

X! 4) 

as 

X4 

X! 5 ) 

s 

X5 

Y ( 1 ) 

= 

Y 1 

Y( 2) 

= 

Y2 

Y< 3 ) 

= 

Y3 

Y ( 4) 

= 

Y4 

Y( 5 ) 

= 

Y5 

DO 11 


K = 1,5 

IF! Y( K ) )1 2, 1 1 , 1 1 

Y(K ) 


AB S ( Y ( K ) ) 

NAG = 


1 


I 1 CONTINUE 

IF(NAG) 13,60,13 
13 WRITE! 6,6 ) 

6 FORMA T ( 77HJNEGA TI VE Y VALUES HAVE BEEN MADE POSITIVE IN CALCUL AT 13 
INS BEYOND THIS LINE.. ///> 

60 DO 130 1=1,5 
H{ I ) - ALOG! Y( I ) ) 

130 CONTINUE 

ZUCK = O.ODO 
ZACK = O.ODO 
7.0 RP = O.ODO 
AKE = O.ODO 
AKK = O.ODO 
ARP = O.ODO 
DO 1C I = 1,5 
ZUCK = ZUCK + W ( I ) 

ZACK = ZACK +■ X! I )*W< I) 

ZORP = ZORP + H ( I )*w( I ) 

AKE = AKE + X( I )*W(I) 

AKK = AK< + W{ I )*X< I )**2 
10 ARP = ARP + X( I )*H( I )*W( I ) 

FARGO = ZUCK*AKK - ZACK*AKE 
IF! FARGO) 30,20,30 
2^ WR I TE ( 6, 8 ) 

8 FORMA T(30H DIVISION BY ZERO NO SOLUTION.) 

R E T UR N 

30 APPLE = ( ZORP * AKK - Z AC K* AR P) ZF A RG 0 
ALPHA =DEXP( APPLE) 

BETA = ( Z UCK*ARP - ZORP*AKE ) /FARGO 

1000 Y IN T = ALPHA * EXP (B: TA*X! NT) 

RET UR N 
END 



$ IEF iC SCRE 


DECK 


SUBROUTINE SPIN I X , Y » N , Z , YI NT , MA X > 

C SPLINE ROUTINF 

DIMENSION X( 256 ) ,YI2S6J ,2(256) , YINTI 2 56) 

DOUBL E PRECISION S(256),AI256)e8I256),CI256},F!256),Wf256),SBI256! 
1 1 G ( 256),EM(256),0Y0X<256) ,D2YDX(25fe) ,CURV<255 ) ,RA0I 256)»SJM(?56) 

DO 10 I = 2? N 

10 S( I ) = X( I ) - X(I-l) 

NO = N - 1 

I F( 2.GT.N0 ) GO TO 25 

00 20 I = 2, NO 

A( n = S( l ) /6.0 

B(I) = ( S( I )+S( IHII /3.0 

c< i ) = si 1+1 > /6.o 

20 F(I) = (Y(!tl)-Y(III/SIIH)-(YIII-Y(I-lll/Sm 
2 5 AIN) = -0.50 
B< 1 )= 1.0 
B I N > = 1.0 
Cl 1 ) =-0.50 
Fill =0.0 
FI N 1 = 0.0 
Will = B( 1) 

sbi l) = c ( n/wm 

G( 1 ) = 0.0 

00 30 I = 2»N ^ 

W( I) = BI Il-A (I )* SB 1 1- L I 
sbi 1 1 = c m / wi i > 

30 G< I ) = IF(lH(ll*Gll-in/W([| 

EMIN) = GIN) 

DO AO I = 2,N 
K = N +1-1 

AO EMIK) = GIK) - SB I K ) *E M I K+ 1 ) 

SUM! 1 )= 0.0 
DO 45 I = 2,N 

4 5 SUM! I ) =SUM I I- 1I+SI l>* I Y(I>+Yt 1-1 > ) /2.-SII ) **3*( EMI I ) + FMt I- 1 ) )/24. 
DO 90 I = 1,MAX 
K = 2 

IF ( ZI I >-X ( 1) ) 60,50,70 
50 YINTI I )=Y( 1) 

GO TO 66 

60 IF I Z I I ) .L T . X I 1 ) ) YINTI I ) = 0.0 
IF! Z( I I.L T.XI 1) ) GO TO 90 
GO TO 8 5 

65 IF! ZI I J.GT.XIN) ) Y I NT 1 1 ) = 0.0 
IF! ZI I J.GT.XIN) ) GO TO 90 

K = N 
GO TO 6 5 

70 I F ( ZI I )- X I K ) ) 85, 75,80 
75 YINTI I ! = YI K 1 
GO TO 66 
80 K = K +1 

IF! K-N ) 70, 70,65 

8 5 Y INTI I )=EM(K- 1 )* I X I K ) - Z 1 1 )) **3/6. /S(K) +EMIKI *( Z I I ) -XIK-l ) ) **3/ 6 . 
l/SI K ) +( YI K )/S(K)-FM(K)*S<K) /6. ) * (ZI I ) -X(K-l)f + lY(K-l)ZSIK) - EMI K-l ) 
2* SI K ) / 6. )*( XI K )-Z II ) ! 

86 OYOXI I)=-EM(K-l)*IX<K)-Z(m**2/2.0/SIK) + EM(K>*lX(K-l)-ZII))**2/2. 
IO/SIK )+IY (K)-Y(K-l) ) /S(K)-(EM(K)— EMIK-l) ) *S(K)/6.0 
0 2YDXI U = EM!K-1)*(XIK5-ZII ) )/S(K)+EM(K)*(Z( I ) -X ( K- 1 ) ) / S C K) 
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Figure 1 . - Portion of output display of two-dimensional analyzer. 


(Ub 4W41) 
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Relative neutron flux, number/(cirr)(sec)(MeV) 



Figure 2. - Series of pulse height distributions obtained using stilbene 
and polonium-210 - beryllium source (ref. 2). 



Neutron energy, MeV 

Figure 3. - Plutonium beryllium spectra unfolded by FERDOR input prepared manually and by PREJUD. 


30 


NASA- Langley, mi 


22 


E-5824 


National Aeronautics and Space Administration 
Washington, D. C. 20546 


official business 


FIRST CLASS MAIL 


POSTAGE AND FEES PAID 
NATIONAL AERONAUTICS AND 
SPACE ADMINISTRATION 


POSTMASTER: 


If Undeliverable (Section 158 
Postal Manual) Do Not Return 


'’The aeronautical and space activities of the United States shall be 
conducted so as to contribute ... to the expansion of human knowl- 
edge of phenomena in the atmosphere and space. The Administration 
shall provide for the widest practicable and appropriate dissemination 
of information concerning its activities and the results thereof." 

— National Aeronautics and Space Act of 1958 


NASA SCIENTIFIC AND TECHNICAL PUBLICATIONS 


TECHNICAL REPORTS: Scientific and 
technical information considered important, 
complete, and a lasting contribution to existing 
knowledge. 

TECHNICAL NOTES: Information less broad 
in scope but nevertheless of importance as a 
contribution to existing knowledge. 

TECHNICAL MEMORANDUMS: 
Information receiving limited distribution 
because of preliminary data, security classifica- 
tion, or other reasons. 

CONTRACTOR REPORTS: Scientific and 
technical information generated under a NASA 
contract or grant and considered an important 
contribution to existing knowledge. 


TECHNICAL TRANSLATIONS: Information 
published in a foreign language considered 
to merit NASA distribution in English. 

SPECIAL PUBLICATIONS: Information 
derived from or of value to NASA activities. 
Publications include conference proceedings, 
monographs, data compilations, handbooks, 
sourcebooks, and special bibliographies. 

TECHNOLOGY UTILIZATION 
PUBLICATIONS: Information on technology 
used by NASA that may be of particular 
interest in commercial and other non-aerospace 
applications. Publications include Tech Briefs^ 
Technology Utilization Reports and 
Technology Surveys. 


Details on the availability ol these publications may be obtained from: 

SCIENTIFIC AND TECHNICAL INFORMATION OFFICE 

NATIONAL AERONAUTICS AND SPACE ADMINISTRATION 

Washington, D.C. 20546 


